home *** CD-ROM | disk | FTP | other *** search
/ The Fatted Calf / The Fatted Calf.iso / Applications / Audio / Spectro / Source / SignalProcessor.m < prev    next >
Text File  |  1992-01-22  |  10KB  |  428 lines
  1.  
  2. /* Generated by Interface Builder */
  3.  
  4. #import "SignalProcessor.h"
  5. #import <stdio.h>
  6. #import <math.h>
  7. #import <string.h>
  8. #define SRATE 22050
  9. #define PI 3.141592654782
  10. #define SQRT_TWO 1.414213562
  11. #define TWO_PI 6.283185309564
  12.  
  13. float sines[16384],cosines[16384];
  14. int last_length = 0;
  15.  
  16. @implementation SignalProcessor
  17.  
  18. - window: (int) size array: (float *) array    //   call with defaults Hanning and size/2 centered
  19. {
  20.     return [self window: size array: array type: "Hanning" phase: 0];
  21. }
  22.  
  23. - window: (int) size array: (float *) array type: (char *) type phase: (BOOL) phase
  24. //    Window array of floats with window of type <"Triangle","Hanning","Hamming",
  25. //        "Blackman3","Blackman4","Kaiser"> 
  26. //        with window centered at zero (phase=TRUE) or size/2 (phase=FALSE)
  27. {
  28.     extern triangle(int size,float* array,BOOL phase),
  29.     hanning(int size,float* array,BOOL phase),
  30.     hamming(int size,float* array,BOOL phase),
  31.     blackman3(int size,float* array,BOOL phase),
  32.     blackman4(int size,float* array,BOOL phase),
  33.     kaiser(int size,float* array,BOOL phase);
  34.  
  35.     if (!strcmp(type,"Triangle"))         triangle(size,array,phase);
  36.     else if (!strcmp(type,"Hanning"))         hanning(size,array,phase);
  37.     else if (!strcmp(type,"Hamming"))         hamming(size,array,phase);
  38.     else if (!strcmp(type,"Blackman3"))     blackman3(size,array,phase);
  39.     else if (!strcmp(type,"Blackman4"))     blackman4(size,array,phase);
  40.     else if (!strcmp(type,"Kaiser"))         kaiser(size,array,phase);
  41. // add your own windows here!!!
  42.     else fprintf(stdout,"Windowing Rectangular\n");
  43.     return self;
  44. }
  45.  
  46. void triangle(int size,float* array,BOOL phase)
  47. {
  48.     int i;
  49.     float w,win_frq;
  50.     win_frq = TWO_PI / size;
  51.     if (phase)
  52.     for (i=0;i<size/2;i++)    {
  53.         w = size/2 - i;
  54.         array[i] *= w;
  55.         array[size-i] *= w;
  56.     }
  57.     else
  58.     for (i=0;i<size/2;i++)    {
  59.         w = i;
  60.         array[i] *= w;
  61.         array[size-i-1] *= w;
  62.     }
  63. }
  64.  
  65. void hanning(int size,float* array,BOOL phase)
  66. {
  67.     int i;
  68.     float w,win_frq;
  69.     win_frq = TWO_PI / size;
  70.     if (phase)
  71.     for (i=0;i<size/2;i++)    {
  72.         w = 0.5 + 0.5 * cos(i * win_frq);
  73.         array[i] *= w;
  74.         array[size-i] *= w;
  75.     }
  76.     else
  77.     for (i=0;i<size/2;i++)    {
  78.         w = 0.5 - 0.5 * cos(i * win_frq);
  79.         array[i] *= w;
  80.         array[size-i] *= w;
  81.     }
  82. }
  83.  
  84. void hamming(int size,float* array,BOOL phase)
  85. {
  86.     int i;
  87.     float w,win_frq;
  88.     win_frq = TWO_PI / size;
  89.     if (phase)
  90.     for (i=0;i<size/2;i++)    {
  91.         w = 0.54 + 0.46 * cos(i * win_frq);
  92.         array[i] *= w;
  93.         array[size-i] *= w;
  94.     }
  95.     else
  96.     for (i=0;i<size/2;i++)    {
  97.         w = 0.54 - 0.46 * cos(i * win_frq);
  98.         array[i] *= w;
  99.         array[size-i] *= w;
  100.     }
  101. }
  102.  
  103. void blackman3(int size,float* array,BOOL phase)
  104. {
  105.     int i;
  106.     float w,win_frq;
  107.     win_frq = TWO_PI / size;
  108.     if (phase)
  109.     for (i=0;i<size/2;i++)    {
  110.         w = (0.42 + 0.5 * cos(i * win_frq) + 0.08 * cos(2 * i * win_frq));
  111.         array[i] *= w;
  112.         array[size-i] *= w;
  113.     }
  114.     else
  115.     for (i=0;i<size/2;i++)    {
  116.         w = (0.42 - 0.5 * cos(i * win_frq) + 0.08 * cos(2 * i * win_frq));
  117.         array[size-i] *= w;
  118.         array[i] *= w;
  119.     }
  120. }
  121.  
  122. void blackman4(int size,float* array,BOOL phase)
  123. {
  124.     int i;
  125.     float w,win_frq;
  126.     win_frq = TWO_PI / size;
  127.     if (phase)
  128.     for (i=0;i<size/2;i++)    {
  129.         w = (0.35875 + 0.48829 * cos(i * win_frq) + 
  130.             0.14128 * cos(2 * i * win_frq) + 0.01168 * cos(3 * i * win_frq));
  131.         array[i] *= w;
  132.         array[size-i] *= w;
  133.     }
  134.     else
  135.     for (i=0;i<size/2;i++)    {
  136.         w = (0.35875 - 0.48829 * cos(i * win_frq) + 
  137.             0.14128 * cos(2 * i * win_frq) - 0.01168 * cos(3 * i * win_frq));
  138.         array[i] *= w;
  139.         array[size-i] *= w;
  140.     }
  141. }
  142.  
  143. void kaiser(int size,float* array,BOOL phase)
  144. {
  145.     blackman4(size,array,phase);            //  I'm Cheating, see Harris paper, eq. 34
  146. }
  147.  
  148. - logMag: (int) size array: (float *) f
  149. {
  150.     [self logMag: size array: f floor: -60.0];
  151.     return self;
  152. }
  153.  
  154. - log: (int) size array: (float *) f floor: (float) fl
  155. {
  156.     int i;
  157.     float t1=0.0,t=0.0,t2;
  158.     t2 = fl / 10.0;
  159.     if (f[0] > t1)
  160.         t1 = f[0];
  161.     for (i=1;i<size;i++)    {
  162.     if (f[i]>t1)
  163.        t1 = f[i];
  164.     }
  165.     for (i=0;i<size;i++)    {
  166.         t = log10(f[i]/t1);
  167.     if (t<t2)
  168.        t = t2;
  169.     f[i] = 1 - t/t2;
  170.     }
  171.     return self;
  172. }
  173.  
  174. - log: (int) size array: (float *) f floor: (float) fl ceiling: (float) ceiling
  175. {
  176.     int i;
  177.     float t1,t=0.0,t2;
  178.     t2 = fl / 10.0;
  179.     t1 = ceiling;
  180.     for (i=0;i<size;i++)    {
  181.         t = log10(f[i]/t1);
  182.     if (t<t2)
  183.        t = t2;
  184.     f[i] = 1 - t/t2;
  185.     }
  186.     return self;
  187. }
  188.  
  189. - logMag: (int) size array: (float *) f floor: (float) fl
  190. {
  191.     int i;
  192.     float t1=0.0,t=0.0,t2=-12;
  193.     t2 = fl / 10.0;
  194.     f[0] = 2 * fabs(f[0]);
  195.     if (f[0] > t1)
  196.         t1 = f[0];
  197.     for (i=1;i<size/2;i++)    {
  198.         t = f[i]*f[i] + f[size-i]*f[size-i];
  199.     f[i] = t;
  200.     if (t>t1)
  201.        t1 = t;
  202.     }
  203.     printf("%f\n",t1);
  204.     for (i=0;i<size/2;i++)    {
  205.         t = log10(f[i]/t1);
  206.     if (t<t2)
  207.        t = t2;
  208.     f[i] = 1 - t/t2;
  209.     }
  210.     return self;
  211. }
  212.  
  213. - logMag: (int) size array: (float *) f floor: (float) fl ceiling: (float) ceiling
  214. {
  215.     int i;
  216.     float t1=0.0,t=0.0,t2=-12;
  217.     t2 = fl / 10.0;
  218.     f[0] = 2 * fabs(f[0]);
  219.     for (i=1;i<size/2;i++)    {
  220.         t = f[i]*f[i] + f[size-i]*f[size-i];
  221.     f[i] = t;
  222.     }
  223.     t1 = ceiling;
  224.     for (i=0;i<size/2;i++)    {
  225.         t = log10(f[i]/t1);
  226.     if (t<t2)
  227.        t = t2;
  228.     f[i] = 1 - t/t2;
  229.     }
  230.     return self;
  231. }
  232.  
  233. - magnitude: (int) size array: (float *) array magOut: (float *) mag
  234. {
  235.     int i;
  236.     mag[0] = 2 * fabs(array[0]);
  237.     for (i=1;i<size/2;i++)    {
  238.         mag[i] = sqrt(array[i] * array[i] + array[size-i] * array[size-i]);
  239.     }
  240.     return self;
  241. }
  242.  
  243. - square_mag: (int) size magnitudeIn: (float *) mag_in squareMagOut: (float *) square_mag_out
  244. {
  245.     int i;
  246.     for (i=0;i<size;i++)  square_mag_out[i] = mag_in[i] * mag_in[i];
  247.     return self;
  248. }
  249.  
  250. - normalize: (int) size array: (float *) array;
  251. {
  252.     int i;
  253.     float mx = 0.0;
  254.     for (i=0;i<size;i++) if (fabs(array[i])>mx) mx=fabs(array[i]);
  255.     if (mx!=0.0) for (i=0;i<size;i++) array[i] /= mx;
  256.     return self;    
  257. }
  258.  
  259. - dht: (float *) input_array output: (float *) output_array length: (int) length
  260. {
  261.     int i,j;
  262.     float w;
  263.     w = TWO_PI / length;
  264.     for (i=0;i<length;i++)    {
  265.         output_array[i] = 0.0;
  266.     for (j=0;j<length;j++)    {
  267.         output_array[i] += input_array[j] * (cos(w * i * j) + sin(w * i * j));
  268.     }
  269.     }
  270.     return self;
  271. }
  272.  
  273. void make_sines(int length)
  274. {
  275.     int i;
  276.     float freq,temp;
  277.     if (length!=last_length)    {
  278.         last_length = length;
  279.         freq = 2.0 * PI / length;
  280.         for (i=0;i<length;i++) {
  281.             temp = freq * i;
  282.         sines[i] = sin(temp);
  283.         cosines[i] = cos(temp);
  284.         }
  285.     }
  286. }
  287.  
  288. -  fhtRX2: (int) powerOfTwo array: (float *) array
  289. {
  290.     return self;        //   Not implemented yet
  291. }
  292.  
  293. - fhtRX4: (int) powerOfFour array: (float *) array
  294. {
  295.     /*  In place Fast Hartley Transform of floating point data in array.
  296.         Size of data array must be power of four. Lots of sets of four 
  297.     inline code statements, so it is verbose and repetitive, but fast. 
  298.     A 1024 point FHT takes approximately 80 milliseconds on the NeXT computer
  299.     (not in the DSP 56001, just in compiled C as shown here).
  300.  
  301.     The Fast Hartley Transform algorithm is patented, and is documented
  302.     in the book "The Hartley Transform", by Ronald N. Bracewell.
  303.     This routine was converted to C from a BASIC routine in the above book,
  304.     that routine Copyright 1985, The Board of Trustees of Stanford University     */
  305.     
  306.     #define PI 3.141592654782
  307.     #define SQRT_TWO 1.414213562
  308.     #define TWO_PI 6.283185309564
  309.  
  310.     register int j=0,i=0,k=0,L=0;
  311.     int n=0,n4=0,d1=0,d2=0,d3=0,d4=0,d5=1,d6=0,d7=0,d8=0,d9=0;
  312.     int L1=0,L2=0,L3=0,L4=0,L5=0,L6=0,L7=0,L8=0;
  313.     int n_over_d3;
  314.     float r=0.0;
  315.     int a1=0,a2=0,a3=0;
  316.     float t=0.0,t1=0.0,t2 =0.0,t3=0.0,t4=0.0,t5=0.0,t6=0.0,t7=0.0,t8=0.0;
  317.     float t9=0.0,t0=0.0;
  318.     n = pow(4.0 , (double) powerOfFour);
  319.     if (n!=last_length)    {
  320.         make_sines(n);
  321.     last_length = n;
  322.     }
  323.     n4 = n / 4;
  324.     r = SQRT_TWO;
  325.     j = 1;
  326.     i = 0;
  327.     while (i<n-1)    {
  328.         i++;
  329.     if (i<j)    {
  330.             t = array[j-1];
  331.         array[j-1] = array[i-1];
  332.         array[i-1] = t;
  333.         }
  334.         k = n4;
  335.         while ((3*k)<j)    {
  336.             j -= 3 * k;
  337.         k /= 4;
  338.         }
  339.         j += k;
  340.     }
  341.     for (i=0;i<n;i += 4) {
  342.         t5 = array[i];
  343.     t6 = array[i+1];
  344.         t7 = array[i+2];
  345.     t8 = array[i+3];
  346.     t1 = t5 + t6;
  347.         t2 = t5 - t6;
  348.         t3 = t7 + t8;
  349.         t4 = t7 - t8;
  350.     array[i] = t1 + t3;
  351.     array[i+1] = t1 - t3;
  352.     array[i+2] = t2 + t4;
  353.     array[i+3] = t2 - t4;
  354.     }
  355.     for (L=2;L<=powerOfFour;L++)  {
  356.     d1 = pow(2.0 , L+L-3.0);
  357.     d2=d1+d1;
  358.     d3=d2+d2;
  359.     n_over_d3 = n / 2 / d3;
  360.     d4=d2+d3;
  361.     d5=d3+d3;
  362.     for (j=0;j<n;j += d5)      {
  363.          t5 = array[j];
  364.          t6 = array[j+d2];
  365.          t7 = array[j+d3];
  366.          t8 = array[j+d4];
  367.          t1 = t5+t6;
  368.          t2 = t5-t6;
  369.          t3 = t7+t8;
  370.          t4 = t7-t8;
  371.          array[j] = t1 + t3;
  372.          array[j+d2] = t1 - t3;
  373.          array[j+d3] = t2 + t4;
  374.          array[j+d4] = t2 - t4;
  375.          d6 = j+d1;
  376.          d7 = j+d1+d2;
  377.          d8 = j+d1+d3;
  378.          d9 = j+d1+d4;
  379.          t1 = array[d6];
  380.          t2 = array[d7] * r;
  381.          t3 = array[d8];
  382.          t4 = array[d9] * r;
  383.          array[d6] = t1 + t2 + t3;
  384.          array[d7] = t1 - t3 + t4;
  385.          array[d8] = t1 - t2 + t3;
  386.          array[d9] = t1 - t3 - t4;
  387.          for (k=1;k<d1;k++)    {
  388.           L1 = j + k;
  389.           L2 = L1 + d2;
  390.           L3 = L1 + d3;
  391.           L4 = L1 + d4;
  392.           L5 = j + d2 - k;
  393.           L6 = L5 + d2;
  394.           L7 = L5 + d3;
  395.           L8 = L5 + d4;
  396.           a1 = (int) (k * n_over_d3) % n;
  397.           a2 = (a1 + a1) % n;
  398.           a3 = (a1 + a2) % n;
  399.           t5 = array[L2] * cosines[a1] + array[L6] * sines[a1];
  400.           t6 = array[L3] * cosines[a2] + array[L7] * sines[a2];
  401.           t7 = array[L4] * cosines[a3] + array[L8] * sines[a3];
  402.           t8 = array[L6] * cosines[a1] - array[L2] * sines[a1];
  403.           t9 = array[L7] * cosines[a2] - array[L3] * sines[a2];
  404.           t0 = array[L8] * cosines[a3] - array[L4] * sines[a3];
  405.           t1 = array[L5] - t9;
  406.           t2 = array[L5] + t9;
  407.           t3 = - t8 - t0;
  408.           t4 = t5 - t7;
  409.           array[L5] = t1 + t4;
  410.           array[L6] = t2 + t3;
  411.           array[L7] = t1 - t4;
  412.           array[L8] = t2 - t3;
  413.           t1 = array[L1] + t6;
  414.           t2 = array[L1] - t6;
  415.           t3 = t8 - t0;
  416.           t4 = t5 + t7;
  417.           array[L1] = t1 + t4;
  418.           array[L2] = t2 + t3;
  419.           array[L3] = t1 - t4;
  420.           array[L4] = t2 - t3;
  421.          }
  422.     }
  423.     }          
  424.     return self;
  425. }
  426.  
  427. @end
  428.